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Electronic  Structure  Methods  based  on  Density  Functional  Theory 

Christopher  Woodward 

Materials  and  Manufacturing  Directorate 

Air  Force  Research  Faboratory 

Introduction 

Over  the  last  two  decades  electronic  structure  methods,  based  on  Density  Functional  Theory, 
have  emerged  as  a  powerful  tool  for  assessing  the  mechanical,  thermodynamic  and  defect 
properties  of  metal  alloys.  These  “First  Principles”  methods  are  very  appealing  because  they  are 
based  on  the  culmination  of  our  understanding  of  quantum  mechanics  and  the  electron-ion  many- 
body  problem.  While  the  starting  point  for  such  calculations  requires  only  the  most  a  basic 
knowledge  of  chemistry,  crystalline  and  defect  structure  the  calculations  can  quickly  become 
very  computational  challenging  with  increasing  system  size  and  complexity.  Practical 
application  of  electronic  structure  methods  invariably  includes  chemical,  spatial  or  temporal 
approximations  that  can  curtail  a  faithful  representation  of  the  actual  materials  problem. 
However,  over  the  last  decade  there  have  also  been  significant  advances  in  methods  for 
calculation  free  energies  (entropy)[l],  activated  states  (i.e.  kinetics)[2],  flexible  boundary 
conditions [3],  lattice  dynamics[4],  and  reaction  rate  theory[5].  Taken  with  the  rapid 
improvements  in  computer  processor  speeds  and  the  maturation  of  “easy-to-use”  DFT  methods 
there  has  been  an  explosive  growth  in  the  use  of  DFT  methods  in  materials  science. 

This  Chapter  is  meant  to  be  a  guide  to  understanding  the  origins  of  these  methods,  their 
strengths  and  limitations  and  provide  the  basic  procedures  for  calculating  essential  structural 
properties  in  metal  alloys.  Before  delving  into  the  details  for  DFT  it  is  important  to  place  the 
method  in  the  context  of  the  larger  community  of  scientists  interested  in  the  nature  of  the 
electronic  state  in  materials.  Modern  electronic  structure  emerged  from  method  development  in 
the  Chemistry  and  the  Physics  communities  and  fall  roughly  into  two  groups:  Hartree  Fock  and 
its  extensions  (HF+E)  and  Density  Functional  Theory.  Historically,  HF+E  was  been  considered 
to  be  more  precise,  and  has  been  preferred  by  Physical  Chemists,  because  systematic 
improvements  to  the  original  approximation  are  well  defined.  For  several  reasons  HF+E  is  not 
well  suited  for  metallic  systems,  as  will  be  discussed  in  detail  in  section  N.l,  and  corrections  to 
HF  are  extremely  computationally  intensive,  scaling  with  the  4th  to  6th  power  of  the  number  of 
electrons.  Density  Functional  Theory  has  been  widely  used  in  metallic  systems  since  it’s 
inception  in  1962.  With  improvements  in  efficiency  (speed)  and  refinements  in  the  underlying 
approximations  DFT  is  increasingly  being  used  in  Quantum  Chemistry  applications.  Recently 
researchers  have  begun  to  blur  the  line  between  these  two  approaches  by  constructing  novel 
potentials  blending  fundamental  aspects  of  the  two  theories  [6].  The  resulting  approximations 
show  great  promise  for  calculations  over  a  broad  range  of  problems  ranging  from  atoms  and 
molecules  to  chemically  complex  metal-oxide  interfaces. 

For  scientists  and  engineers  considering  using  electronic  structure  methods,  navigating  the  sea  of 
DFT  acronyms  can  be  challenging.  In  general  the  acronyms  refer  to  the  numerical  scheme,  or 
basis,  used  to  represent  the  electrons.  More  recently,  as  methods  have  matured  codes  have  been 
named  after  the  groups  that  developed  or  support  the  method.  All  electronic  structure  methods 
must  deal  with  the  large  changes  in  the  electron  distribution  observed  in  atoms,  molecules  and 
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solids.  Some  of  the  electrons  are  strongly  bound  to  the  nuclear  sites  (core  states)  and  are  very 
similar  to  that  found  in  isolated  atoms.  Others  are  more  weakly  bound  (valence  states)  producing 
bonding,  and  are  responsible  for  most  of  the  electronic,  optical,  thermodynamic  and  chemical 
properties.  Electronic  structure  methods  deal  with  this  disparity  in  a  variety  of  ways,  depending 
on  what  class  of  materials  problem  is  under  consideration.  For  example  with  isolated  atoms  and 
molecules  it  is  natural  to  work  in  real  space  with  methods  based  on  a  linear  combination  of 
atomic  orbitals  (LCAO),  represented  numerically  or  as  a  sum  of  analytic  functions  (e.g. 
gaussians).  For  crystalline  systems  it  is  more  natural  to  use  periodic  boundary  conditions  and 
electrons  are  represented  using  a  linear  combination,  or  basis,  of  plane  waves.  Over  time  several 
methods  were  developed  to  avoid  the  large  number  of  planewaves  needed  to  represent  the 
rapidly  varying  electron  core  densities.  One  approach,  employed  in  augmented  plane  wave 
(APW)  and  muffin  tin  orbital  (MTO)  methods,  is  to  use  a  set  of  local  functions  centered  around 
each  atom  and  to  match  that  solution  on  a  sphere  to  a  plane  wave  basis  everywhere  else.  Another 
technique,  the  pseudopotential  method,  maps  the  strongly  bound  electron  states  into  a  potential 
that  is  then  used  to  calculate  the  valence  electrons.  These  Pseudopotential  Plane  Wave  (PPW) 
methods  are  relatively  easy  to  use  and  the  simplicity  of  the  basis  has  allowed  significant  progress 
in  computational  efficiency  (i.e.  parallel  processing)  and  analytic  solutions  for  properties  such  as 
atomic  forces  and  stress. 

Historically  mixed  basis  methods  (APW,  MTO)  have  been  considered  as  the  benchmark  for 
accuracy  in  most  applications.  However,  the  mixed  basis  makes  these  methods  more  challenging 
to  use  and  adds  significant  complications  to  deriving  basic  quantities  such  as  atomic  forces  or 
stress.  With  advances  in  pseudopotential  theory  since  the  mid  1980’s  PPW  methods  routinely 
reproduce  the  results  of  mixed  basis  methods.  Also,  because  of  the  added  benefits  of  larger 
simulation  sizes,  automated  atomic  and  cell  optimization  it  can  be  argued  PPW  methods  are 
producing  more  accurate  results  is  wider  range  of  applications. 

Mixed  basis  methods  are  still  the  preferred  technique  for  systems  where  the  pseudopotential 
approximation  breaks  down.  This  happens  when  changes  in  valence  electrons  change  the 
structure  of  the  core  electrons.  For  example  in  the  actinides,  where  the  f-electrons  are  coupled  to 
the  core  states,  APW  and  MTO  methods  are  preferred[7].  Also,  simulations  of  photo  absorption 
and  emission  are  probably  best  modeled  using  techniques  where  the  core  states  are  optimized 
along  with  the  valence  states.  The  quality  and  availability  of  pseudopotentials  in  some  software 
packages  can  be  quite  limited.  Researchers  new  to  the  field  should  carefully  assess  the  available 
options  before  investing  time  or  resources  in  any  particular  method. 

One  additional  criterion  to  consider  is  the  scale  of  the  problem  that  needs  to  be  solved.  Inevitably 
this  dictates  the  method  and  the  required  computational  platform.  Both  FCAO  and  PW 
pseudopotential  methods  scale  well  on  current  parallel  supercomputers  and  are  typically  applied 
to  molecular  and  crystalline  problems  respectively. 

DFT  is  being  applied  throughout  the  scientific  community  to  a  staggering  range  of  problems. 
With  current  multi-processor  supercomputers  PPW  methods  can  simulate  system  sizes  up  to 
approximately  1000  atoms[8].  This  varies  with  the  system  symmetry  and  the  choice  of  atomic 
species,  with  the  transition  metals  being  the  most  challenging.  Researchers  are  also  running  ab- 
initio  molecular  dynamics  simulations  for  cells  ranging  from  100-500  atoms  for  simulation  times 


2 


up  to  tens  of  pico-seconds.  The  future  of  DFT  will  be  driven  by  improvements  to  the  underlying 
approximations,  the  introduction  of  new  hybrid  potentials,  and  advances  in  supplementary 
methods  the  use  employ  DFT  results.  Research  into  new  novel  basis  functions  (e.g.  wavelets),  or 
the  introduction  of  new  computer  hardware  (e.g.  Field  Programmable  Gate  Arrays)  could 
revolutionize  the  field.  Finally,  while  still  in  it’s  infancy  there  is  a  significant  effort  underway  to 
directly  calculate  the  electronic  state  by  Quantum  Monte-Carlo  methods,  if  properly  coupled  to 
next  generation  supercomputers  this  eventually  could  overtake  all  other  developments  [9]. 

The  rest  of  this  chapter  is  divided  into  three  independent  sections.  The  first  reviews  the  general 
underlying  theory  of  electronic  structure  methods  and  Density  Functional  Theory  specifically 
and  the  taxonomy  of  DFT  methods  that  have  emerged  over  the  last  thirty  years.  The  second 
section  reviews  the  approximations  and  computational  details  of  the  most  popular  method  used 
in  metal  systems,  the  pseudopotential  plane  wave  methods.  The  last  section  reviews  a  subset  of 
the  applications  of  DFT  methods  have  found  in  metals  alloy  systems.  This  includes  calculations 
of  a  variety  of  structural,  thermodynamic  and  defect  properties  with  particular  emphasis  on 
structural  metal  alloys  and  their  derivatives. 

N.l  The  Fundamentals  of  Density  Functional  Theory: 

We  would  like  to  model  a  chunk  of  matter  using  only  what  we  know  about  coulomb  interactions 
between  electrons  and  ions  and  the  underlying  principles  of  quantum  mechanics.  The  approach 
taken  over  the  last  50  years  has  been  to  systematically  apply  approximations  making  the  many- 
body  problem  more  manageable  while  retaining  the  essential  physics.  Part  of  this  evolving 
approach  is  to  reduce  the  systems  of  equations  to  that  subset  which  captures  the  problem  of 
interest.  We  are  not  interested  in  solving  systems  with  Avogadro’s  number  of  particles;  not  only 
would  solving  such  a  problem  be  unfeasible,  analyzing  the  results  of  such  a  calculation  would  be 
a  herculean  task.  Therefore  for  practical  reasons,  both  conceptual  and  computational,  it  is 
considered  best  practice  to  minimize  the  scale  (spatial  and  temporal)  of  the  electronic  structure 
calculation.  There  are  many  good  reviews  Density  Functional  methods.  For  a  general  overview 
of  the  fundamentals  see  R.  Martin’s  text,  Electronic  Structure,  Basic  Theory  and  Practical 
Methods [10].  Payne  et  al.  review  the  general  theory  of  Pseudopotential  Plane  Wave  methods  and 
many  practical  issues  on  applying  these  methods [1 1].  Many  practical  details  of  PPW  and 
Augmented  Planewave  methods  are  reviewed  in  D.  Singh’s  text  Planewaves,  Pseudopotentials 
and  the  Linearized  Augmented  Plane  Wave  Method[12]. 

Beginning  from  classical  mechanics,  the  many  body  Hamiltonian  of  an  ensemble  of  interacting 
atoms  takes  the  form: 
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Where  Mi,  Pi,  Zi  and  Ri  are  the  mass,  momentum,  charge  and  position  of  the  M  possible  ions  and 
me,  pe  and  r,  are  the  mass,  momentum  and  position  of  the  m  possible  electrons.  The  Hamiltonian 
is  then  separated  into  two  parts,  a  purely  ionic  part  (the  first  two  terms  in  equation  1)  and  an  ion- 
electron  part: 

Hrotai  ~  HIon(.Ri)  +  He_Ion(RI,ri)  (2) 

From  this  point  Materials  Scientist  can  choose  to  represent  He-ion  by  an  effective  potential  which 
leads  to  the  field  of  atomistic  modeling,  or  to  invoke  Quantum  Mechanics  and  solve  the  many 
body  Schrodinger  equation  which  leads  to  the  field  of  electronic  structure  methods.  The  later  are 
generally  considered  a  more  faithful  representation  of  the  many  body  problem,  as  the  electrons 
are  treated  explicitly.  Standard  practice  is  to  now  apply  the  Bom-Oppenheimer  approximation. 
Given  that  the  electron  cloud  responds  much  faster  to  an  applied  field  than  the  ions  (me/Mi«l), 
we  can  decouple  the  nuclear  and  electronic  motion  and  solve  for  the  electron  degrees  of  freedom 
with  the  ionic  positions  held  fixed.  Using  separation  of  variables  the  Schrodinger  equation 
corresponding  to  equation  1  can  be  divided  into  two  parts: 

rj)  =  (He  +  Hlonm{R„  rj)  =  E  r J)  (3) 

using  =  ^Fe ({/?/,  ri})vF/on ({/?/})  produces  two  equations: 

(Te  +  Ve-e  +  Fe_i)Te( {R„ri})  =  £,evTe({R/,ri})  and  (T,  +  V,_,  +  FJ'F,  ({«,})  = 

£’/'P/0„({R/})  where  T  and  V  refer  to  the  kinetic  and  coulomb  potential  terms  for  the  electrons 
and  ions  and  the  eigenvalues  (Ee  coming  from  the  separation  of  variables)  incorporated  as  an 
effective  potential  for  the  ionic  problem.  The  many  body  Schrodinger  equation  for  the  m 
electrons  is: 


in 

He'¥e({RI,ri})=Yj 


i=  1 


=  EeVe{{RI,ri}) 


({«/,  rj) 


Unfortunately  this  equation  cannot  be  solved  directly.  Two  approaches  are  taken  to  solve  this 
system  of  equations  the  first,  Hartree-Fock  and  it’s  extensions  (HF+E),  solves  for  the  electron 
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wave-functions  and  the  second,  Density  Functional  Theory,  solves  for  the  charge  density  (Drs. 
W.  Kohn  and  J.A.  Pople  split  the  Nobel  Prize  in  Chemistry  in  1998  for  aspects  of  these 
contributions).  The  Hartree-Fock  approach  is  attractive  because  the  derivation  allows  for  well- 
defined  systematic  (though  costly)  improvements  to  the  initial  approximation  for  the  third  term 
in  equation  (2) ,  these  are  sometimes  called  “post  Hartree  Fock  methods”.  HF+E  methods  are 
used  extensively  in  non-metallic  systems,  however  they  are  poorly  suited  for  metal  systems  for 
several  reasons.  First,  significant  corrections  to  the  initial  HF  approximation  are  required  to 
properly  represent  metal  systems.  Second,  in  free  electron  metals  HF  produces  an  intrinsic 
instability  in  the  electron  velocity  (a  logarithmic  divergence  in  ds/dk,  where  e(k)  is  the  energy 
dependence  of  the  electron  as  a  function  of  wave  vector  k)  at  the  Fermi  surface[13]. 

Density  Functional  theory  is  based  on  two  insights  provided  Hohenberg,  Kohn  and  Sham  in  the 
early  1960’s[14,15].  First,  Hohenberg  and  Kohn  proved  that  for  the  ground  state  of  an  interacting 
electron  gas  in  an  external  potential  the  electron  density,  p(r),  can  be  treated  as  the  total  energy 
of  a  system  of  interacting  electrons  in  an  external  potential  (i.e.  the  coulomb  potential  produced 
by  the  atomic  nuclei)  and  is  given  exactly  as  a  functional  of  the  ground  state  electronic  density: 

E  —  E (p) .  Here  a  functional  is  defined  as  a  function  of  a  function  -  in  this  case  E  is  function  of 
the  electron  density.  While  the  Hohenberg-Kohn  theorem  shows  that  E(p)  is  a  unique  functional 
it  does  not  provide  a  prescription  on  how  to  form  the  functional,  so  the  usefulness  of  the  theorem 
is  dependent  on  finding  sufficiently  accurate  approximations  [14]. 

To  this  end  Kohn  and  Sham  suggested  writing  E  as: 

E  [p]  =  Ts  [p]  +  Eei  [p]  +  Eu  [p]  +  Eh  [p]  +  Exc  [p]  (4) 


the  functionals  on  the  right  representing  the  kinetic  energy  of  a  system  of  non-interacting 
electrons,  the  electron-ion  interactions,  the  Hartree  potential  of  electron-electron  interactions,  the 
ion-ion  interactions  and  the  exchange  correlation  functional  respectively[15].  The  functional  are 
now  integrals  over  space,  i.e.  the  last  two  terms  explicitly  are: 


EH[P ]  =  7//  dr  dr 


pO)p(r') 


(5) 


Exc  [p]  =  /  dr  p(r)  6xc(p(r)) 


(6) 


In  the  last  term  Kohn  and  Sham  identified  exc(p(r))  as  the  exchange  and  correlation 
energy/electron  of  a  uniform  electron  gas  of  density  p.  This  is  the  local  density  approximation 
(LDA)  which  assumes  that  given  a  sufficiently  slowly  varying  density  a  function,  exc,  can  be 
defined  which  represents  the  effective  potential  of  an  electron  surrounded  by  its  own  "mutual 
exclusion  zone"  consistent  with  the  requirements  of  quantum  mechanics. 
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Using  the  fact  that  the  functional  E[p ]  is  an  energy  minimum  with  respect  to  variations  in  p  (the 
H-K  theorem)  they  then  derived  single  particle  Schrodinger  like  equations  that  are  sometimes 
referred  to  as  the  Kohn  Sham  equations: 

(  ft?  v~>  1  r  p(v'')  ) 

[2^V2  ~  Ze2  lu\r  -  R,\  +  e2  J  \r  —  r'\  +  ^c(p(r))|  W  =  U0;O)  (7) 

Where  p(r)  =  SI^il0i(r)  I2  with  m  equal  to  the  number  of  occupied  states  (the  number  of 
electrons  in  the  system)  and  pxc(p)  =  d(pexc(p))/dp  which  is  identified  as  the  exchange- 
correlation  contribution  to  the  chemical  potential  of  a  uniform  gas  of  densityp.  The  system  of 
equations  is  solved  self-consistently  by  assuming  a  p(r)  constructing  the  last  two  terms  in  the 
KS  equations  and  then  solving  for  p(r)  using  0,(r).  The  total  energy  is  given  by: 

m 

E[p(r. )]  =  ^  G  -  P^^'\  drdr  +  /  p(r^£xc(p(r^  ~  ftxc{p(r))\dr  (8) 

t= i 

The  K-S  equation  maps  the  interacting  many  electron  system  to  a  set  of  non-interacting  electrons 
moving  in  an  effective  potential  of  all  the  other  electrons.  The  utility  of  the  K-S  equations  rest  in 
our  ability  to  find  reasonable  approximations  for  the  functional  Exc  [p] .  Fortunately  this  function 
has  been  studied  in  detail  for  the  case  of  a  uniform  electron  gas  [16],  and  derived  using  Monte 
Carlo  techniques [17]  and  parameterized  for  electronic  structure  calculations [18]. 

LDA  has  been  surprisingly  successful  in  predicting  a  variety  of  properties  in  metals.  Lattice 
parameters  are  usually  accurate  to  within  ~1%  and  cohesive  energies  and  elastic  constants  to 
within  -10%.  The  method  is  well  suited  for  studying  solids,  perfect  and  defected  crystals  and  is 
easily  extended  to  include  spin  dependence,  the  Local  Spin  Density  Approximation,  which  has 
been  widely  applied  to  ferromagnetic  and  anti-ferromagnetic  systems  [19].  LDA  is  also  the 
starting  point  for  a  variety  of  improvements  based  on  the  local  change  in  the  electron  density 
produced  by  the  electron  (the  exchange-correlation  hole).  These  Generalized  Gradient 
Approximations  have  systematically  improved  the  accuracy  of  DFT  for  problems  in  molecular 
systems  and  broadened  the  application  base  significantly  (these  are  reviewed  in  the  next  section). 
However,  currently  the  method  is  still  not  well  suited  for  systems  with  large  Van  der  Waals 
energies  or  systems  sampling  infinitesimally  small  electron  densities  such  as  a  structures 
bounded  by  a  vacuum. 

Where  DFT  methods  seem  to  diverge  is  in  the  spatial  representation  of  the  one-electron  wave 
functions.  The  wide  variety  of  methods  reflects  the  fact  that  an  accurate  representation  of  the 
charge  density  has  traditionally  required  specialized  basis  functions.  The  fundamental  problem  is 
that  as  the  atomic  number  increases  the  additional  atomic  wave-functions  are  required,  by  the 
Pauli  exclusion  principle,  to  be  orthogonal  to  existing  lower  lying  wave-functions.  To 
accomplish  this,  as  the  principle  quantum  number  increases,  the  wave-functions  take  on  a  rapidly 
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varying  radial  form  near  the  atomic  center.  Therefore,  for  a  set  of  basis  functions  to  accurately 
describe  the  electrons  it  must  be  able  to  both  represent  the  rapidly  varying  function  near  atomic 
centers  and  the  relatively  smooth  functional  form  outside  that  region. 

Several  strategies  have  been  used  to  solve  this  problem.  Techniques  such  as  the  Augmented 
Plane  Wave,  Muffin  Tin  Orbital  and  Korringa  Kohn  Rostoker  methods  use  an  efficient  and 
compact  basis  to  describe  the  wave-functions  near  the  atomic  centers  [20,2 1,22].  An  additional 
basis  (typically  plane-waves)  is  used  to  describe  the  wave-functions  outside  this  region  and 
various  schemes  are  used  to  ensure  a  proper  match  of  the  wave-functions  at  the  boundary 
between  these  two  regions.  These  methods  give  a  very  accurate  representation  of  the  core  region, 
while  allowing  some  flexibility  in  the  basis  outside  this  region.  Refinements  to  these  techniques 
such  as  the  Full  potential  Linearized  Augmented  Plane  Wave  (FLAPW)  [23]  and  Full  Potential 
Linearized  Muffin  Tin  (FP-LMTO)  [24]  are  currently  considered  to  be  the  most  accurate  DFT 
methods,  though  the  implementations  are  limited  to  relatively  small  cell  sizes. 

The  seminal  work  on  Orthogonal  Plane -Wave  methods  led  to  the  development  of  an  alternative 
technique  where  the  low  lying  (core)  states  are  effectively  removed  from  the  calculation [25].  In 
this  case  an  effective  electron-ion  potential,  or  pseudo-potential,  is  derived  from  an  atomic 
electronic  structure  calculation.  Pseudo-potentials  incorporate  the  tightly  bound  wave-functions 
and  ionic  charge  so  that  the  potential  produces  the  same  electronic  interactions  as  the  original 
atomic  calculation.  In  this  way  the  core  electrons,  which  do  not  normally  influence  materials 
properties,  can  be  removed  from  the  simulation.  By  construction  the  pseudo-potentials  produce 
the  same  interaction  with  the  valence  electrons  as  the  original  all  electron  calculation  (as 
measured  through  the  electron  scattering  properties) [26].  Recent  pseudo-potential  schemes  have 
relaxed  this  philosophy  during  the  construction  of  the  potential,  only  to  re-impose  the 
requirement  when  constructing  the  Kohn-Sham  orbitals  in  the  system  of  interest[27,28].  Pseudo¬ 
potentials  also  incorporate  the  effective  potential  produced  by  the  Pauli  exclusion  principle  such 
that  the  valence  wave-functions  are  smooth  functions  in  all  space.  Therefore  it  is  natural  to 
combine  a  plane-wave  basis  with  the  pseudo-potential  representation  in  what  we  now  call 
pseudo-potential  plane- wave  methods  (PPW). 

While  the  plane-wave  basis  is  not  as  compact  as  that  used  in  the  LAPW  methods,  PPW  methods 
have  been  easier  to  implement  because  of  the  simplicity  of  the  plane-wave  basis.  It  is  relatively 
straightforward  to  calculate  atomic  forces  (through  the  Hellmann-Feynman  theorem),  the  stress 
tensor  and  phonon  properties  [29,30,31].  However,  until  the  mid  1980’s  application  of  PPW 
methods  were  somewhat  limited  because  of  the  size  of  the  required  plane-wave  basis  set.  In  1985 
Car  and  Parrinello  showed  how  to  simultaneously  optimize  both  the  electronic  and  ionic  degrees 
of  freedom  by  taking  advantage  of  fast  fourier  transforms  and  the  plane-wave  basis[32].  Iterative 
diagonalization  methods  that  have  grown  out  of  these  insights  have  shown  that  by  proper 
preconditioning  of  the  Kohn-Sham  wave-functions  during  the  optimization  procedure  it  is 
possible  to  directly  minimize  the  Kohn-Sham  energy  functional[33].  This  innovation  stabilizes 
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the  constrained  optimization  of  Eq.  1  and  has  made  it  possible  to  run  very  large  simulations 
relatively  efficiently[ll]. 

Figure  1.  Schematic  of  the  methods  developed  by  the  chemistry  and  physics  communities  to 
solve  the  materials  many  body  problem.  In  general  the  DFT  methods  are  the  method  of  choice 
for  calculations  on  metallic  systems. 
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N.2  Pertinent  Approximations  and  Computational  Detail  for  Calculations  in  Metal  Alloys: 

In  order  to  calculate  the  electronic  structure  of  an  alloy  the  researcher  needs  more  than  just  an 
underlying  theory  and  a  working  description  (basis)  of  the  electrons.  First  the  methods  need 
efficient  and  accurate  techniques  for  integrating  the  quantities  described  in  section  1  over  the 
volumes  (and  k-space)  of  interest.  Second,  the  approximations  and  basis  need  to  be  well  matched 
to  the  problem  of  interest.  Finally,  explicit  knowledge  of  the  crystal  structure,  location  and 
species  of  every  atom  on  the  simulation  volume  is  required.  Fortunately,  for  a  most  of  the  alloy 
of  engineering  interest  this  information  is  available  in  tabulations  of  the  International  Tables  for 
Crystallography  in  the  form  of  space  groups,  describing  the  symmetry  of  the  lattice,  and 
Wyckoff  positions,  describing  the  atomic  sites  (position  and  chemistry) [34]. 
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Integration  of  cell  quantities:  One  advantage  to  working  on  metallic  systems  is  that  the 
underlying  crystalline  structure  is  almost  invariable  periodic.  This  allows  the  researcher  to 
employ  simulation  cells  with  periodic  boundary  conditions  (supercells)  to  represent  the  material 
of  interest.  Periodic  boundary  conditions  also  make  it  possible  to  represent  most  all  quantities  of 
interest  in  terms  of  real  and  Fourier  space  (k- space)  components.  This  is  particularly  useful  when 
summing  up  terms  numerically  in  the  KS  equations  or  equation  (8)  over  the  Brillouin  zone. 

The  most  straightforward  way  to  integrate  of  quantities  over  the  supercell  is  to  divide  the 
Brillouin  zone  using  a  tetrahedron  grid  of  points  (k-points).  However,  in  the  early  1970’s 
numerical  schemes  were  developed  to  predict  the  smallest  set  of  k  points,  “special  k-point”,  that 
would  yield  the  most  accurate  cell  integrations  [35,36].  In  most  modern  application  codes  the 
selection  of  k-points  is  sufficiently  automated  that  the  investigator  needs  to  only  input  the 
required  density  of  such  points.  However,  one  characteristic  of  a  metal  is  that  valence  states,  or 
bands,  are  partially  occupied,  and  this  can  produce  numerical  instabilities.  In  order  to  avoid  this 
and  to  improve  the  efficiency  of  the  cell  integration  researchers  introduced  a  numerical  smearing 
of  the  highest  lying  bands,  specifically  those  that  are  near  the  Fermi  surface.  The  “broadening” 
methods  fall  into  two  classes,  those  employing  ad-hoc  functional  forms  such  as  Gaussians  and 
finite  temperature  schemes  based  on  the  Fermi-Dirac  or  Gauss  like  functions  that  mimic  thermal 
broadening  [37,38].  Both  methods  are  effective,  and  the  latter  technique  has  the  advantage  of  an 
associated,  if  ad-hoc,  temperature. 

Computational  time  for  metal  simulations  scales  with  the  number  of  k-points,  however  as  the 
simulation  sizes  get  larger  and  reciprocal  space  get  smaller  the  number  of  required  k-points  is 
reduced.  For  large  simulations,  say  greater  than  100  atoms,  often  only  a  single  k-point  is  needed. 
Also,  if  only  the  gamma  point  (k  =  (0,0,0))  is  required  then  calculations  can  gain  another  factor 
of  2  in  efficiency,  because  of  the  symmetry  imposed  on  the  complex  parts  of  the  wavefunctions. 

While  there  are  rules  of  thumb  for  the  use  of  special  k-points  and  broadening  methods,  it  pays  to 
carefully  test  the  convergence  of  such  methods  using,  for  example,  the  cell  energy  or  other 
quantity  of  interest.  In  general  the  density  of  k-points  should  be  approximately  that  of  the  size  in 
dispersion  in  the  bands  near  the  Fermi  surface[39]. 

Understanding  and  choosing  a  pseudopotential:  Early  in  the  development  of  electronic 
structure  methods  researchers  realized  that  the  electrons  contained  in  full  atomic  shells  (s,  p  and 
d)  do  not  have  a  strong  influence  on  chemical  and  mechanical  properties.  These  effects  are 
controlled  mostly  by  the  interaction  of  valence  electrons,  which  have  the  largest  principle 
quantum  number,  and  thus  the  most  rapidly  varying  radial  function  in  the  region  around  the  atom 
nucleus.  Pseudopotentials  replace  the  core-valence  electron  interactions,  the  second  term  of 
equations  (4)  and  (7),  with  an  effective  potential  produces  a  realistic  pseudo-valence 
wavefunction  that  has  a  smooth  and  slowly  varying  radial  form.  Typically  pseudopotentials  are 
derived  from  an  atomic  reference  calculation  and  then  used  in  crystalline  or  other  environment, 
so  transferability  is  a  serious  concern  when  developing  such  a  scheme.  Early  local 
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pseudopotentials,  and  more  recent  implementations  of  the  same,  are  severely  limited  and  can 
only  reliable  be  used  in  simple  free  electron  metals  where  the  core  electrons  have  very  weak 
interactions  with  the  valence  states[40.41]. 

Modern  pseudopotential  theory  is  based  on  a  “norm-conservation”  approach  that  enforces  a  strict 
criterion  for  mapping  real  to  pseudo  wavefunctions  and  includes  non-local  angular  momentum 
(£)  dependent  interactions  that  accurately  model  the  valence-core  electron  interactions  [26]. 
Transferability  is  maintained  by  imposing  identical  logarithmic  derivatives,  and  thus  scattering 
phase  shifts,  outside  a  certain  (core)  radius  about  each  atomic  site.  More  recent  advances  in 
pseudopotential  theory  such  as  Vanderbilt’s  ultrasoft  pseudopotentials  make  use  of  additional 
functions  about  the  atomic  core,  which  allow  for  smoother  pseudo  wavefunctions  and  more 
efficient  PPW  calculations [27].  This  and  later  refinements  of  pseudopotential  theory  parallel  the 
original  strategies  used  in  Orthogonalized  Plane  Wave  methods  developed  by  Herring,  Callaway 
and  others  from  the  mid  1900’s[25].  The  most  recent  advances  in  pseudopotential  theory,  the 
projector  augmented  wave  (PAW)  method[28],  retains  all  the  information  of  the  core  states  and 
is  thus  analogous  to  the  most  accurate  all  electron  methods  (e.g.  OPW,  APW  and  MTO). 
Implementing  the  PAW  methods  in  PPW  codes  required  additional  development  and  using  the 
potentials  incurs  additional  computational  overhead. 

Using  modern  numerical  methods,  current  commercial  PPW  implementations  of  the  PAW 
method  are  as  efficient  as  the  original  Car  and  Parrinello  methods  and  as  accurate  as  many  full 
potential  methods.  They  have  the  added  benefit  of  ease  of  calculation  of  atomic  forces,  stress 
tensor  and  convergence  of  basis.  A  wide  variety  of  ultrasoft  and  PAW  pseudopotentials  are 
available  in  the  user  community  as  well  as  source  code  for  developing  such  potentials.  More 
importantly  there  are  readily-available,  well-documented  suites  of  pseudopotentials  that  have 
been  tested  by  a  broad  user  base. 

Exchange- Correlation  Potentials,  Local  Density  Approximation  and  the  Generalized 
Gradient  Approximation:  Density  Functional  Theory  is  one  of  the  most  successful  electronic 
structure  methods  precisely  because  of  the  simplicity  of  the  underlying  exchange-correlation 
functional.  Practical  application  of  the  HK  theorem  through  the  KS  equations  requires  a  both  an 
assessment  of  the  exchange  correlation  functional  and  a  numerically  efficient  scheme  for 
interpolating  the  energy  for  a  range  of  charge  densities.  Since  the  original  formulation  of  the  KS 
equations  the  nature  of  the  exchange  correlation  potential  has  been  studied  in  some  detail.  The 
Local  Density  Approximation  is  the  foundation  of  all  these  approximations.  Within  the  LDA 
only  knowledge  of  the  exchange  correlation  energy  of  the  homogeneous  electron  gas  is  required. 
Thiis  approximated  as  the  sum  of  exchange  and  correlation  potentials,  the  first  given  by  a  basic 
analytical  form  and  the  second  calculated  using  Monte  Carlo  methods  [17].  These  data  were  then 
fit  to  functional  forms  to  improve  computational  efficiency  and  parameterized  for  electronic 
structure  calculations  [18]. 
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The  Local  Density  Approximation  has  been  found  to  be  a  surprisingly  accurate  in  a  wide  variety 
of  systems.  The  initial  formulation  was  expected  only  to  valid  for  volumes  with  slowly  varying 
electron  densities,  a  condition  that  is  not  well  satisfied  in  many  crystals.  It  is  generally  believed 
that  the  LDA  approximation  underestimates  the  exchange  energy  (by- 10%)  and  overestimates 
the  correlations  energy  (2%)  and  that  these  errors  partially  cancel  each  other  out  [42].  However, 
LDA  is  particularly  unsatisfactory  for  low  electron  densities,  such  near  a  surface,  and  that  has 
made  the  approximation  problematic  for  calculations  of  atoms  and  molecules.  Still,  LDA 
produces  reasonable  accurate  bond  lengths  and  geometries  for  some  molecules. 

The  efficacy  of  the  KS  equations  and  the  need  for  highly  accurate  simulations  has  resulted  in 
systemic  improvements  to  the  LDA.  The  most  successful  approaches,  based  on  generalized 
gradient  approximations  (GGA),  include  information  on  the  effects  of  inhomogeneities  in  the 
electron  gas  on  the  exchange  correlation  potential.  The  gradient  corrections  are  constructed  to 
satisfy  intrinsic  sum  rules  and  are  designed  to  maintain  the  accuracy  of  LDA  while  correcting  the 
errors  introduced  by  large  gradients.  Using  FP-LMTO  Ozolins  and  Korling  calculated  the 
changes  in  lattice  constants  and  bulk  modulus  produced  by  using  the  GGA  proposed  by  Perdew 
and  Zhang  (sometimes  referred  to  as  PW91)[43].  They  found  a  systematic  improvement  in 
equilibrium  volumes  and  bulk  modulus  for  3d,  4d,  and  5d  transition  metals,  with  the  mean  error 
decreasing  on  average  by  50%  for  both  quantities  across  the  series[44].  Other  researchers  also 
found  that  early  GGA  methods  and  PW91  correctly  produces  the  correct  bee  ground  state  for 
crystalline  Fe  where  the  Local  Spin  Density  Approximation  erroneously  predicts  an  fee  ground 
state  [45.46]. 

Recently,  other  GGA  methods  are  being  validated  that  produce  better  energetics  and  better 
represent  low  density  regions.  Hybrid  schemes  based  on  a  weighted  mixing  of  HF  exchange  and 
DFT  correlation  have  gained  favor  in  the  quantum  chemistry  community  [47].  Also,  another 
class  of  GGA  functional  has  been  self  consistently  matched  to  high  and  low  electron  densities, 
making  it  efficient  and  well  suited  for  metallic  systems  with  internal  or  external  surfaces  [48] 

N.3  Practical  Application  of  DFT  in  Metals  and  Alloys: 

As  illustrated  in  Figure  N.2  DFT  methods  come  in  a  variety  of  forms,  they  also  vary  widely  in 
their  level  of  maturity  and  efficacy.  Until  recently  most  of  the  mixed  basis  methods  (FPLMTO 
and  FLAPW)  were  closely  held  academic  codes,  now  there  are  several  freeware  and  commercial 
options  (see  Table  3).  While  highly  accurate,  these  methods  have  a  more  complex  set  of 
adjustable  basis  parameters  than  PPW  methods.  There  are  a  variety  of  PPW  methods  some 
available  as  freeware  and  some  commercially,  and  some  have  well  developed  pseudopotential 
libraries.  New  users  should  verify  a  given  code  base  has  the  necessary  features  for  running  their 
application  before  investing  resources  into  a  method.  All  of  the  methods  shown  in  Table  3  should 
produce  accurate  results  for  the  simple  applications  outlined  in  this  section. 
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Table  3.  Partial  list  of  currently  available  DFT  programs.  Some  are  available  for  a  license  fee  and 
others  are  available  at  no  cost.  Many  of  the  methods  have  associated  users  groups  and  some  have 
graphical  user  interfaces.  Materials  Science  trade  organizations  are  beginning  to  track  the  status 
of  software  (see  for  example:  http://iweb.tms.org/forum)  and  may  provide  useful  updates  to  this 
table. 


Method 

Acronym/Name 

Fee 

POC 

ABINIT 

No 

httn://www. abinit.org/.  Prof.  X.  Gonze, 
Universite  catholique  de  Louvain,  Physico- 
Chemistry  and  Physics  of  Materials,  Louvain- 
la-Neuve,  BELGIUM 

PPW 

CASTEP 

Yes 

Accelrys,  Inc.,  San  Diego,  CA  92121 
httD://accelrvs.com/Droducts/materials-studio 

Quantum  Expresso:  opEn 
Source  Package  for  Research 
in  Electronic  Structure, 
Simulation,  and  Optimization 

No 

P.  Giannozzi,  Universit'a  di  Udine  and 
Democritos  National  Simulation  Center,  Italy 
http  ://www  .quantum-espresso  .org 

VASP:  Vienna  Ab-initio 
Simulation  Package 

Yes 

Prof.  J.  Hafner,  Institute  of  Materialphysik 

Wien  University  Austria, 

http  ://cms .  mpi.  univie .  ac .  at/vasp/ 

LCAO 

DM0L3 

Yes 

Accelrys,  Inc,  San  Diego,  CA  92121 
http://accelrvs.com/products/materials-studio 

SIESTA 

Yes 

Prof.  J.A.  Torres,  Universidad  Autonoma  de 

Madrid,  Spain 

http://www.icmab.es/siesta, 

http://www.nanotec.es/ 

FP- 

LMTO 

LmtART 

No 

http://www.fkf.mpg.de/andersen/ 

S.  Y.  Savrasov,  Phys.  Rev.  B  54,  16470  (1996). 

RSPt :  Relativistic  Spin 
Polarised  (test) 

No 

http :  //w  ww  .rspt  .net/index .  php 

WIEN2K 

Yes 

http:// www. wien2k. at/inde x . html ,  Prof. 

Karlheinz  Schwarz,  Inst.  f.  Materials 

Chemistry,  TU  Vienna 

FP- 

LAPW 

FLAIR 

Prof.  M.  Wienert,  Univ.  Wisconsin  Milwaukee, 
weinert@uwm.edu 

QMD-FLAPW 

Yes 

Prof.  A.J.  Freeman,  Northwestern  Univ. 

Quantum  Materials  Design,  Inc., 
http :  //flap  w .  com/ne  w  s  .html 

From  crystal  structure  to  input  file,  examples  of  VASP  input  files: 

Setting  up  the  simulation  cell  for  electronic  structure  calculations  requires:  the  lattice  vectors,  the 
atomic  positions  and  the  chemical  species  at  each  site.  The  input  file  required  for  the  PPW 
method  VASP  will  be  used  to  illustrate  this  process  [49,50].  Typically  the  researcher  starts  with 
a  phase  and  then  refers  to  tables  and  textbooks  to  find  the  required  quantities.  For  example  in  Ni- 


12 


based  superalloys  Ni3Al  is  an  important  precipitate  that  significantly  strengthens  the  Ni  matrix 
phase.  Using  the  tabulated  data  on  alloy  phases  in  W.E.  Pearson’s  Handbook  of  Lattice  Spacings 
and  Structures  of  Metals  [51]  the  structure  type  is  listed  having  a  structure  name  (a  representative 
material)  AuCu3,  the  Strukturbericht  designation  LI2,  a  lattice  parameter  of  3.567  Angstrom, 
with  a  space  group  of  Pm3m.  The  Pearson  classification  for  LI2,  AuCu3  is  given  as  cP4  in  the 
tables  leading  up  to  the  Table  of  Classification  of  Structures  of  Metals  and  Alloys.  A  shortened 
version  of  the  entry  under  cP4  is  given  in  Table  4.  After  the  first  line  describing  the  structure 
designations  the  atomic  basis,  the  Wyckoff  positions,  are  listed  in  terms  of  the  atom  type,  the 
number  of  atoms  at  each  symmetry  distinct  point  and  the  internal  coordinate.  All  the  required 
information  is  now  determined.  The  International  Tables  for  X-Ray  crystallography  has  more 
information  for  this  cubic  space  group  (number  221)  and  lists  much  more  complicated  crystal 
structures  with  this  space  group[34,52]. 

Table  4.  Crystallographic  information  and  atomic  basis  for  LI 2,  NI3AI. 


Classification 

symbol 

Structure  name 

Strukturbericht  „ 

type  SpaCC  gr°Up 

cP4 

AuCu3 

LI  2  Pm3m 

Origin  at  center  (m3m) 
Equivalent  positions: 
Au:  1  a  m3m 

0,0,0 

Cu:  3  c  4/mmm  0,  W/2;  Vi, 0,  W,  VM). 

The  VASP  the  input  file  for  the  lattice  vectors  and  atomic  basis  is  called  POSCAR  and  a  screen 
shot  of  the  POSCAR  file  for  Ni3Al  is  shown  below.  In  POSCAR  the  title  line  is  followed  by  the 
lattice  constant  from  Pearson,  the  three  lattice  vectors  for  a  cubic  lattice  in  Cartesian  coordinates. 
This  is  followed  by  the  number  of  each  chemical  species,  in  this  case  3  Ni  atoms  and  one  A1 
atom,  and  a  keyword  describing  the  format  of  the  atomic  basis.  (Note  that  the  order  of  the 
chemical  species  is  important  and  must  be  consistent  with  the  ordering  of  the  pseudopotential 
input  file.)  The  atomic  basis  can  be  entered  in  either  Cartesian  (Keyword:  CARTESIAN)  or  in 
terms  of  the  three  lattice  vectors  (Keyword:  DIRECT).  For  the  “Direct”  mode  the  atomic 
positions  correspond  to  R  =  x1a1  +  x2a2  +  x3 a3  ,  where  aL  are  the  lattice  vectors  scaled  by  the 
lattice  parameter  and  Xi  are  the  values  entered  into  POSCAR.  Input  using  the  CARTESIAN 

— )  /s 

keyword  is  scaled  only  by  the  lattice  constant:  R  =  xti  +  x2j  +  x3k  ,  where  i,  j,  and  k  are  unit 
vectors  [100],  [010],  and  [001]  respectively. 

Fig.  2.  Screen  shot  of  input  file  describing  crystal  system  and  atomic  basis  for  y’-NLAl. 
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Example:  Ni3Al 
3,57640 

1,00000  0.00000  0.00000 
0.00000  1.00000  0.00000 
0.00000  0.00000  1.00000 
3  1 

Direct 

0.000  0.500  0.500 
0.500  0.000  0.500 
0.500  0.500  0.000 
0.000  0.000  0.000 


Title 

lattice  Constant  (Ang.) 

Lattice  vector 
Lattice  vector 
Lattice  vector 

Number  of  each  atomic  species 
Format  of  the  atomic  positions 
Internal  coordinates  of  atoms 


Another,  more  complex  example  of  an  atomic  basis  is  the  5-MoNi  phase,  which  is  also  important 
in  the  Ni-based  superalloys.  In  Pearson  this  topologically  closed  packed  phase  is  listed  as 
orthorhombic  with  lattice  constants:  a=9.108,  b=9.108,  and  c=  8.852  Ang.,  containing  56  atoms 
with  the  space  group  P2 1 2 1 2 1 .  Table  4  gives  the  representative  atomic  positions  listed  as  14 
roman  numbers  with  corresponding  coordinates.  The  Wyckoff  positions,  listed  in  the 
International  tables,  for  P2i2i2j,  are  shown  in  Table  5.  The  atomic  coordinates  are  generated  by 
using  the  last  column  of  Table  5  with  each  of  the  14  atomic  parameters  producing  the  expected 
56  atoms. 

Table  4.  Crystallographic  information  and  atomic  parameters  for  5-MoNi  [51]. 


Phase 

System 

Strukturbericht 

type 

Space 

group 

Est.  % 
Mo[53] 

5-MoNi 

Ortho¬ 

rhombic 

-None- 

P2i2i2i  Atoms 

Atomic  parameters 

IV 

0.4519 

0.1153 

0.5322 

0 

VI 

0.4424 

0.3662 

0.5972 

0 

VIII 

0.3882 

0.0523 

0.2748 

0 

IX 

0.1337 

0.0707 

0.2157 

0 

X 

0.3768 

0.4358 

0.8567 

0 

XII 

0.0680 

0.1442 

0.9529 

0 

I 

0.1763 

0.4832 

0.6425 

80 

XIII 

0.0338 

0.3398 

0.1807 

80 

II 

0.2289 

0.2865 

0.4098 

91 

V 

0.2648 

0.1993 

0.7486 

91 

XI 

0.3136 

0.2464 

0.0740 

58 

VII 

0.0029 

0.1969 

0.6767 

100 

XIV 

0.1885 

0.0157 

0.4960 

100 

III 

0.1031 

0.4192 

0.9133 

100 

Note  however  the  atomic  species  for  each  site  is  still  unknown,  going  back  to  the  original 
reference  for  this  crystallographic  assessment  we  can  find  the  chemical  assignments  for  most  but 
not  all  the  atomic  sites [53].  At  finite  temperatures  all  alloys  show  deviations  from  perfect 
ordering,  however  the  composition  of  the  studied  5-MoNi  phase  was  Mo49.2-Ni  and  X-Ray 
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analysis  could  not  unambiguously  determine  the  chemistry  on  at  least  one  of  the  sites.  The 
electronic  structure  calculations  can  proceed  by  assuming  a  Mo50-Ni  composition  and  the 
chemistry  at  the  sites  designated  by  XI  as  being  occupied  by  Ni  atoms.  This  is  however  just  one 
possibility,  and  in  principle  a  free  energy  model  would  include  sampling  the  formation  energy  of 
other  atomic  arrangements  at  this  composition. 


Table  5.  Wyckoff  Positions  for  Space  group  P2i2i2i 


Multiplicity 

Wyckoff 

letter 

Site 

Symmetry 

Coordinates 

4 

a 

1 

(x,y,z),(-xU/2,-y,zU/2),(-x,yU/2,-zU/2),(xU/2,-y+V2,-z) 

Using  Tables  4  and  5  the  initial  cell  was  constructed  as  shown  in  the  screen  shot  of  the  POSCAR 
file  in  Figure  3.  The  Figure  also  shows  a  screen  shot  of  the  final  cell  configuration.  Using  VASP 
the  lattice  vectors  and  atomic  positions  were  optimized  within  a  spin  polarized  (ferro-magnetic) 
ultrasoft  pseudopotential  approximation.  Note  the  slight  change  in  the  length  of  the  orthorhombic 
lattice  vectors  and  the  change  in  atomic  positions.  At  the  start  of  the  calculation  the  pressure  was 
~  16  kB  and  the  largest  force/atom  was  ~  0.5  eV/Ang,  after  optimization  the  pressure  was  less 
than  0.5  kB  and  the  atomic  forces  less  than  2e-3  eV/Ang  and  the  total  energy  change  from  the 
intial  configuration  was  ~  0.5  eV. 
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Figure  3.  Screen  shot  of  initial  and  final  cell  configurations  for  5-Mo50-Ni. 


Initial  Configuration  Final  Configuration 


0“  Mo-Ni  (Ni )24  (Ni4  Mol6)  (Mo)12  ** 

9*108 

1*00000  0*00000000  0*00000000 

0*  Mo-Ni  (Ni)24  (Ni4  Mol6)  (Mo)12  ** 

9.10800000000000054 

.9996694028707155  .0000000000000000  .0000000000000000 

0*00000000 

1.00000000 

0.00000000 

*0000000000000000 

*9997687828277412  *0000000000000000 

0*00000000 

0.00000000 

0*97190000 

*0000000000000000 

*0000000000000000  *9767010090926945 

28  28 

Direct 

*45190001 

*11530000 

*53219998 

28  28 

Direct 

*4506643768099344 

,1193068679416401 

*5325479584880161 

*04809999 

*88470000 

*03219998 

*0493356231900656 

,8806931320583600 

*0325479584880161 

*95190001 

*38470000 

*46780002 

*9506643768099343 

,3806931320583599 

*4674520415119839 

*54809999 

*61530000 

*96780002 

.5493356231900657 

,6193068679416400 

*9674520415119839 

.44240001 

*36620000 

*59719998 

*4390199503536315 

,3770221180583645 

*5972503663884575 

.05759999 

*63380000 

*09719998 

*0609800496463684 

,6229778819416355 

*0972503663884576 

*94240001 

*13380000 

.40280002 

*9390199503536316 

,1229778819416355 

*4027496336115424 

*55759999 

*86620000 

*90280002 

*5609800496463684 

,8770221180583645 

*9027496336115425 

*38820001 

*05230000 

*27480000 

.3880875526506487 

,0559391388777245 

.2720464605775263 

*11179999 

*94770000 

*77480000 

*1119124473493513 

,9440608611222755 

*7720464605775261 

*88820001 

*44770000 

*72520000 

*8880875526506486 

,4440608611222755 

.7279535394224739 

*61179999 

*55230000 

*22520000 

.6119124473493514 

,5559391388777245 

*2279535394224737 

*13370000 

*07070000 

*21570000 

*1305762636276348 

,0701291674759629 

*2171473356559112 

*36630000 

*92930000 

*71570000 

*3694237363723653 

,9298708325240371 

*7171473356559112 

*63370000 

*42930000 

*78430000 

*6305762636276347 

,4298708325240371 

*7828526643440888 

.86630000 

*57070000 

*28430000 

*8694237363723653 

,5701291674759629 

*2828526643440889 

*37680000 

*43579999 

*85670000 

*3708537893684525 

,4412869934178103 

*8579881136158205 

*12320000 

*56420001 

*35670000 

*1291462106315474 

,5587130065821898 

*3579881136158205 

*87680000 

*06420001 

.14330000 

*8708537893684526 

,0587130065821897 

*1420118863841867 

*62320000 

*93579999 

*64330000 

*6291462106315474 

,9412869934178102 

*6420118863841795 

*06800000 

*14420000 

*95289999 

*0672333579610334 

,1441663236300990 

*9548821936676196 

*43200000 

*85580000 

*45289999 

*4327666420389666 

,8558336763699008 

*4548821936676197 

*56800000 

*35580000 

*04710001 

.5672333579610335 

,3558336763699009 

*0451178063323803 

*93200000 

*64420000 

*54710001 

.9327666420389665 

,6441663236300992 

.5451178063323804 

*10310000 

*41920000 

*91329998 

*0979147884406859 

,4168979366954007 

*9114546709612028 

*39690000 

*58080000 

*41329998 

*4020852115593141 

,5831020633045993 

*4114546709612028 

*60310000 

*08080000 

*08670002 

*5979147884406859 

,0831020633046064 

*0885453290387972 

*89690000 

*91920000 

*58670002 

*9020852115593141 

,9168979366954007 

*5885453290387972 

*26480001 

*19930001 

*74860001 

.2625803246569123 

,2039912353823866 

.7483085557329274 

*23519999 

.80069999 

*24860001 

*2374196753430877 

,7960087646176135 

*2483085557329273 

.76480001 

*30069999 

*25139999 

*7625803246569123 

,2960087646176135 

*2516914442670727 

*73519999 

*69930001 

.75139999 

*7374196753430877 

,7039912353823865 

*7516914442670726 

*31360000 

*24640000 

*07400000 

.3069190985937987 

,2528840709561812 

*0728398268099706 

*18640000 

*75360000 

*57400000 

*1930809014062013 

,7471159290438190 

*5728398268099705 

.81360000 

*25360000 

*92600000 

*8069190985937985 

,2471159290438189 

*9271601731900295 

.68640000 

*74640000 

*42600000 

*6930809014062015 

,7528840709561810 

*4271601731900294 

*00290000 

*19690000 

*67670000 

*0014759749891245 

,1925737429781877 

*6752655257190752 

*49710000 

*80310000 

*17670000 

*4985240250108754 

,8074262570218123 

*1752655257190751 

*50290000 

*30310000 

*32330000 

*5014759749891246 

,3074262570218124 

*3247344742809249 

*99710000 

*69690000 

*82330000 

*9985240250108754 

,6925737429781877 

*8247344742809248 

*18850000 

*01570000 

*49599999 

*1860527493132290 

,0124308173658511 

*4945098491188892 

*31150000 

*98430000 

*99599999 

*3139472506867710 

.9875691826341489 

.9945098491188892 

*68850000 

*48430000 

*50400001 

*6860527493132291 

,4875691826341489 

.5054901508811108 

.81150000 

*51570000 

*00400001 

.8139472506867709 

,5124308173658511 

*0054901508811108 

*17630000 

*48320001 

.64249998 

*1761110489577723 

,4854782446104694 

*6408034101715941 

*32370000 

*51679999 

*14249998 

*3238889510422276 

,5145217553895307 

*1408034101715943 

*67630000 

*01679999 

*35750002 

*6761110489577723 

,0145217553895306 

*3591965898284058 

*82370000 

*98320001 

*85750002 

.8238889510422277 

,9854782446104693 

.8591965898284059 

*03380000 

*33980000 

*18070000 

*0287512823268599 

,3372034257831562 

.1809562826145831 

*46620000 

*66020000 

*68070000 

*4712487176731401 

,6627965742168439 

*6809562826145832 

*53380000 

*16020000 

*81930000 

*5287512823268599 

,1627965742168439 

*8190437173854168 

*96620000 

*83980000 

*31930000 

*9712487176731401 

,8372034257831561 

.3190437173854170 

*22890000 

*28650001 

*40979999 

*2265797861926567 

,2846050956865158 

*4105987224600145 

*27110000 

*71349999 

*90979999 

*2734202138073433 

,7153949043134844 

*9105987224600145 

*72890000 

*21349999 

*59020001 

*7265797861926566 

,2153949043134842 

*5894012775399855 

*77110000 

*78650001 

*09020001 

*7734202138073434 

,7846050956865156 

*0894012775399856 

This  and  several  other  configurations  were  used  to  develop  a  simple  free  energy  model  of  the  Ni- 
Mo  system  by  approximating  the  configuration  entropy  [54]. 

Lattice  parameters:  While  current  PPW  codes  can  optimize  supercell  geometries  by 
minimizing  in  the  diagonal  components  of  the  stress  tensor,  it  is  still  useful  to  know  how  to 
calculate  the  lattice  parameters  using  equations  of  state.  Perhaps  the  most  cited  equation  of  state 
used  for  this  purpose  was  developed  by  F.D.  Murnaghan  in  1944  [55].  Assuming  that  the  bulk 
modulus  (K)  is  a  linear  function  of  the  pressure: 
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K  =  —V—  =  C(1  +  kP) 


(9) 

Such  that  the  bulk  modulus  (K)  and  its 
derivative  at  zero  pressure 
respectively  are  identified  as: 

C  and  (—) 

\dV  J p= o  dP  V  ctV/ p— o 

Ck  respectively.  Integrating  equation 
(1)  from  zero  pressure  gives  yields: 

P  —  —  —  we  can  identify  K  —  V  — ; 

dV  J  dV2 

and  finally: 


vl£P  =  ^  =  C<!  +  fes((7)“  -  J)  =  C(7)“  <10> 

Integrating  two  times  and  identifying  C  =  K0  and  Ck  =  K0  yields  Mumaghan’s  equation  of  state 
(MES): 


E(V)  =  E(Vo)  + 


VKq 

k'o(Xo- l) 


VqKq 
K'o- 1 


(11) 


This  somewhat  complicated  form  can  then  be  used  to  fit  the  energy  as  a  function  of  volume 
calculated  from  electronic  structure  methods.  For  example  Figure  4  shows  the  energy  vs.  volume 
for  fee  Ni  using  PAW  pseudopotentials  and  Table  6  gives  the  constants  that  produced  the  fitted 
curves. 


In  general  DFT  will  predict  low  temperature  lattice  constants  to  within  a  percent.  The  FDA  will 
typically  underestimate  lattice  parameters,  while  the  later  improved  gradient  corrected 
approximations  do  not  follow  this  trend[57].  MES  also  produces  an  estimate  of  the  bulk 
modulus  and  its  derivative  with  respect  to  volume.  There  are  reasonable  alternatives  to  MES 
such  as  the  Birch  form  [58]  that  is  favored  in  some  applications  [59]. 

Recent  improvements  to  the  gradient  corrections  that  are  designed  to  alleviate  problems  in  low 
charge  density  regions  (i.e.  internal  voids,  surfaces  and  surfaces  interactions)  have  produced 
mean  errors  in  lattice  parameters  of  approximately  0.1%  over  a  wide  range  of  materials  [59]. 
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Table  6.  Results  of  MES  fit  to  VASP  total 

-20.0 

>  -20.5 

>> 

O) 

|  -21.0 

LD 

IB 

o  -21.5 

-22.0 

35  40  45  50  55 

Volume  (A3) 

Figure  4:  Variation  in  fcc-Ni  cell  energy  as  a  Elastic  Constants:  One  of  the  primary  uses  of 

function  of  volume  calculated  using  VASP.  DFT  in  crystalline  metals  has  been  to  predict 

lattice  parameters  and  elastic  constants.  Initially 
the  calculations  were  used  to  assess  the  validity  of  the  LDA  and  the  computational  methods. 
However,  as  the  LDA  became  more  established  and  the  exchange-correlation  functionals  became 
more  refined  it  became  routine  for  groups  to  predict  the  Cij  of  simple  metals.  In  the  early  1990’s 
DFT  was  used  extensively  to  predict  the  elastic  constants  of  a  variety  of  high  temperature 
intermetallics.  Mehl  and  co-workers  at  the  Naval  Research  Laboratory  were  one  of  the  first 
groups  to  apply  these  methods  and  developed  a  robust  strategy  for  assessing  Cij  for  cubic  and 
tetragonal  crystal  structures [60]. 

The  general  approach  is  to  express  the  free  energy  of  the  system  as  a  function  of  the  strain  tensor 
acting  on  a  small  simulation  cell  volume.  We  can  start  from:  dF  —  —SdT  +  PdV  +  dW,  where 
dW  is  the  infinitesimal  work  done  by  elastically  distorting  the  crystal.  Specifically  dW  = 

(Tijdeij  =  Cijki£kid£ij  where  we  have  used  the  definition  of  the  elastic  constants  relating  the 
applied  stress  to  the  resulting  strain:  ct£;-  =  CijkiekL.  Assuming  reversible  and  isothermal  loading 
at  zero  pressure:  dF  =  dW  =  C£/k£6k£d6£/-  and  we  can  write:  d2F / de^de^  =  Cijkl.  Changing 
notations  from  the  fourth  rank  tensor  to  the  reduced  2nd  rank  tensor  [60]  we  express  the  energy  of 
the  system  around  equilibrium  using  a  Taylor  series  expansion  in  the  strain: 

E(e)  =  E0-  P(V0)dV  +  ^f=iZ;6=1  Cyeiej  +  0(e3)  (12) 

V  and  P(V)  are  the  volume  and  pressure  of  the  undistorted  lattice,  dV  is  the  change  in  volume 
produced  by  the  strain  e£.  It  is  natural  to  apply  strains  to  the  simulation  cell  by  transforming  the 
primitive  lattice  vectors  of  the  cell  using  the  strain  tensor  a: 

ai  aT  ei  e6/2  £5/2 

a2  =  a2  (l  +  f)  withs=  66/2  e2  £4/2 

_a3J  LasJ  [es/2  ej2  e3 
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considering  only  non-rotating  strains.  Now  for  the  specific  case  of  cubic  crystals  where: 


C11C12C12  0  0  0 
C12C11C12  0  0  0 
r  r  r  0  0  0 

Cn  —  12  12  11  then  the  double  summation  in  the  equation  for  the  energy  as  a 

1  0  0  0  Q4  0  0  M 

0  0  0  0  C44  0 

0  0  0  0  0  C44- 

function  of  strain  becomes: 

2f=i£y=i  Cije^j  =  ( el  +  ef  +  ef)Cn  +  2(e1e2  +  ^63  +  e2e3  )C12  +  (e|  +  e|  +  e|)C44 

The  effects  of  some  applied  strain  are  now  explicitly  coupled  by  elastic  constants  to  changes  in 
energy  for  a  simulation  cell  with  unit  vectors  ar  For  example,  take  the  case  of  a  hydrostatic 
stress  of  =  e2  =  e3  =  5,  this  yields:  E(e)  =  E0  +  y<52(3C1;L  +  6C12).  Identifying  the  bulk 

modulus  ( K  =  (C41  +  2C12)/3)  we  find:  ^(e)  =  Fq  +  ^V0KS2.  Applying  this  form  to  the  data  in 

Figure  MES_PAW  for  the  spin  polarized  case  yields  bulk  modulus  of  2.01  Mbar,  in  good 
agreement  with  MES. 

In  order  to  find  the  three  independent  elastic 
constants  two  other  equations  are  required,  and 
the  normal  convention  is  to  apply  two  other, 
volume  conserving,  strains.  For  cubic  systems 
usually  C41  —  C12  is  found  by  applying  e±  — 
—e2  =  S,  with  e3  being  set  by  the  constant 

s 2 

volume  constraint  for  a  cubic  cell:  e3  —  — — . 

5  1  -s2 

Similarly,  C44  is  found  by  setting  e 6  =  -  and 

s2 

63  =  4 ^s2'  ^'s  y'clds  two  °lher  equations  for 
the  energy  as  a  function  of  strain: 

E(S)  =  E0  +  E05 2(C14  -  C12)  +  0(54) 
E(S)  =  E0+y52C44+  0(d4) 

The  response  of  a  unit  cubic  cell  of  Ni  to  such 
strains  is  shown  in  Figure  E_CIJ.  When  the  fits  to  the  two  curves  are  combined  with  the  MES 
results  the  elastic  constants  can  be  resolved  as  shown  for  fee  Ni  and  LI2  Ni3Al  in  Tables  2  and  3. 
The  tabulated  results  are  for  spin  averaged,  spin  polarized  (ferro-magnetic),  systems  using  LDA 
and  GGA  approximations.  Tables  7  and  8  show  the  results  from  a  PWPP  calculation  (VASP) 
using  ultrasoft  and  projected  augmented  wave  pseudopotentials  respectively.  Note  that  the  LDA 
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and  LSDA  underestimate  the  lattice  parameter  and  over  estimate  the  elastic  constants  and  the 
GGA  results  (PW91)  show  uniform  improvement  in  lattice  parameters  and  elastic  constants. 
Finally,  though  PAW  is  assumed  to  be  a  better  representation  of  the  core  states  the  USPP’s 
produce  a  more  accurate  misfit  parameter. 

Table  7.  Structural  parameters  for  fee  Ni  and  Ll2  Ni3Al  calculated  using  ultrasoft 
pseudopotentials  in  a  PWPP  method  (VASP).  Spin  averaged  and  spin  polarized  (ferromagnetic) 
calculations  in  the  Local  Density  Approximation  and  a  Generalized  Gradient  Approximation 
(PW91)  are  used  to  predict  the  lattice  parameter  (Angstrom),  elastic  constants  (Mbar)  and  misfit 
parameter  (8  =  2(aNisM  —  aNi) / (aNi  Ai  +  ci-Ni))-  As  expected  the  LDA  and  LSDA 
underestimate  the  lattice  parameters  for  Ni  and  Ni3Al.  GGA  and  SGGA  produce  significantly 
more  precise  lattice  parameters  and  elastic  constants,  with  the  SGGA  calculations  giving  the 
most  accurate  misfit  parameter. 


USPP 

Spin  Averaged 

Spin  Polarized 

Metal 

Property 

LDA 

error 

GGA 

error 

LSDA 

error 

SGGA 

error 

Exp[51 ,561 

Fee 

a0  (A) 

3.4294 

-2.7% 

3.5258 

0.1% 

3.4221 

-2.9% 

3.5337 

0.3% 

3.5238 

Ni 

K 

2.515 

35% 

1.985 

6.8% 

2.383 

28% 

1.962 

5.5% 

1.860 

Cn 

3.154 

27% 

2.506 

1 .0% 

3.034 

22% 

2.380 

-4.1% 

2.481 

C12 

2.195 

42% 

1.725 

11% 

2.057 

33% 

1.753 

13% 

1.549 

C44 

1.358 

9% 

1.057 

-15% 

1.363 

10% 

1.253 

0.9% 

1.242 

<err> 

23% 

6.8% 

19% 

4.8% 

Liz 

a0  (A) 

3.4893 

-2.2% 

3.5769 

0.3% 

3.4928 

-2.1% 

3.5784 

0.3% 

3.5670 

Ni3AI 

K 

2.163 

24% 

1.777 

2.1% 

2.159 

24% 

1.787 

2.7% 

1.740 

Cn 

2.749 

21% 

2.271 

0.3% 

2.778 

23% 

2.397 

5.9% 

2.264 

C12 

1.870 

26% 

1.529 

3.3% 

1.849 

25% 

1.481 

0.1% 

1.480 

C44 

1.431 

11% 

1.173 

-8.7% 

1.471 

15% 

1.240 

-3.4% 

1.284 

<err> 

17% 

2.9% 

18% 

2.5% 

8 

0.0173 

42% 

0.0143 

18% 

0.0204 

68% 

0.0125 

3% 

0.0121 

Table  8.  Structural  parameters  for  fee  Ni  and  Ll2  Ni3Al  calculated  using  Projected  Augmented 
Wave  pseudopotentials  in  a  PWPP  method  (VASP).  Spin  averaged  and  spin  polarized 
(ferromagnetic)  calculations  in  the  Local  Density  Approximation  and  a  Generalized  Gradient 
Approximation  (PW91)  are  used  to  predict  the  lattice  parameter  (Angstrom),  elastic  constants 
(Mbar)  and  misfit  parameter  (8  —  2(aNl.iM  —  aNi) / (aWj3j4i  +  aNi )).  As  expected  the  LDA  and 
LSDA  underestimate  the  lattice  parameters.  For  Ni  and  Ni3Al,  GGA  and  SGGA  produce 
significantly  more  precise  lattice  parameters  and  elastic  constants,  with  the  SGGA  calculations 
giving  the  most  accurate  misfit  parameter. 


PAW 

Spin  Averaged 

Spin  Polarized 

Metal  Property 

LDA  error  GGA  error 

LSDA 

error  SGGA  error 

Exp[51 ,561 
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Fee 

a0  (A) 

3.4197 

-3.0% 

3.5150 

-0.2% 

3.4258 

-2.8% 

3.5219 

-0.1% 

3.5238 

Ni 

K 

1.144 

-38% 

1.968 

5.8% 

1.175 

-37% 

1.942 

4.4% 

1.860 

Cn 

3.185 

28% 

2.480 

0.0% 

3.476 

40% 

2.704 

9.0% 

2.481 

C12 

2.227 

44% 

1.712 

11% 

2.019 

30% 

1.561 

0.8% 

1.549 

C44 

1.383 

11% 

1.121 

-10% 

1.618 

30% 

1.294 

4.2% 

1.242 

<err> 

25% 

5% 

28% 

4% 

L12 

a0  (A) 

3.4823 

-2.4% 

3.5685 

0.04% 

3.4927 

-2.1% 

3.5699 

0.1% 

3.5670 

Ni3AI 

K 

2.183 

25% 

1.779 

2.2% 

2.171 

25% 

1.773 

1 .9% 

1.740 

Cn 

2.774 

23% 

2.264 

0.0% 

2.787 

23% 

2.343 

3.5% 

2.264 

C12 

1.888 

28% 

1.537 

3.9% 

1.863 

26% 

1.488 

0.5% 

1.480 

C44 

1.441 

12% 

1.189 

-7.4% 

1.488 

16% 

1.248 

-2.8% 

1.284 

<err> 

18% 

2.7% 

18% 

1 .8% 

8 

0.0181 

49% 

0.0151 

24% 

0.0193 

59% 

0.0135 

11% 

0.0121 

Entropic  contributions  to  the  free  energy:  In  the  last  10  years  significant  progress  has  been 
made  in  calculating  the  entropic  contributions  to  the  free  energy  of  bulk  phases  and  defects.  This 
includes  configurational,  vibrational  and  electronic  entropic  terms.  Examples  of  applications, 
including  references  reviewing  the  techniques,  are  given  here.  Electronic  entropy  has  been 
shown  to  be  important  in  calculating  defect  energies,  such  as  vacancies  in  body  centered  cubic 
metals  [61].  Contribution  of  thermal  vibrations  to  the  free  energy  as  a  function  of  volume 
(harmonic  and  anharmonic  terms)  has  been  used  to  estimate  the  thermal  expansion  of  a  variety  of 
metals[62,63].  Configurational  entropy  for  dilute  solute  concentrations  are  treated  using  the 
Bragg  Williams  approximations  in  conjunction  with  either  lattice  gas  models  and  the  Low 
Temperature  expansion[64,65].  For  solid  solutions  at  high  concentrations  cluster  expansion 
methods[66,67]  are  used  to  approximate  the  free  energy  on  an  Ising  model  lattice.  Recent 
progress  in  methods  development  has  automated  parts  of  the  construction  and  use  of  these 
techniques[l,68].  Van  de  Walle  and  co-workers  have  also  attempted  to  include  all  three  entropic 
contributions  in  modeling  phases  stability  and  to  inform  CALculation  of  PHAse  Diagram 
methods  (CALPHAD)[69].  These  developments  have  significantly  improved  the  efficiency  and 
accuracy  of  the  Cluster  Expansion  approach,  particularly  in  its  application  to  phase  diagrams. 
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